/* Serial version of the Poisson Equation solver using red-black ordering. */
# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <malloc.h>
# include <math.h>
/*
# include <time.h>
# include <sys/timeb.h>
# include <omp.h>*/

extern void press_bc(const struct grid *g, short tflag, unsigned short iRK);

extern long int solve_poisson(struct grid *g, const double ***d, unsigned short iRK)
  {
   int cnt;
   short i, j, k, nr, nth, nz, im, jm, km, stop_cnt;
   double dpm, dp, ***p, *r, *th, *z, prf;
   char flag;

// struct timeb t1,t2;

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'solve_poisson()'...\n");
#  endif

   r = g->r;   z = g->z;   th = g->th;   nr = g->nr;   nth = g->nth;   nz = g->nz;

   prf = PRF_INIT; stop_cnt = 0;

   p = g->p;

   for(cnt = 1; cnt <= POISS_CNT_MAX; cnt++)
      {
       dpm = 0.0; im = jm = km = -1;

       for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG+(k+i)%2; j < NthG+nth; j += 2)
            {
             dp = P_ijk_NEW - p[i][j][k];

             if((dp<0 ? -dp:dp) > (dpm<0 ? -dpm:dpm))
                  {
                   dpm = dp;
                   im = i; jm = j; km = k;
                  }
             p[i][j][k] += prf*dp;
            }

       for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG+(k+i+1)%2; j < NthG+nth; j += 2)
            {
             dp = P_ijk_NEW - p[i][j][k];

             if((dp<0 ? -dp:dp) > (dpm<0 ? -dpm:dpm))
                  {
                   dpm = dp;
                   im = i; jm = j; km = k;
                  }
             p[i][j][k] += prf*dp;
            }

       press_bc(g, 0, iRK);

       if( (cnt == 1) || (cnt%DISP_FREQ == 0) )
            {
             fprintf(g->dptr, "=> %d: dpm = %E(%d, %d, %d)\n", cnt, dpm, im, jm, km);
             fflush(g->dptr);
            }

       if( ((dpm < 0) ? -dpm:dpm) <= PRESS_DELTA) ++stop_cnt;
       else stop_cnt = 0;

       if(stop_cnt == POISS_CNT_MIN) break;
      }

   flag = (stop_cnt == POISS_CNT_MIN) ? 'Y':'N';

   fprintf(g->dptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c\n", cnt, dpm, im, jm, km, flag);
   fflush(g->dptr);

   fprintf(g->tptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c; ", cnt, dpm, im, jm, km, flag);
   fflush(g->tptr);

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'solve_poisson()' ends.\n");
#  endif

   return cnt;
  }
